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Abstract 

Electrons in one-dimension display the unusual property of separating their spin and charge into 
two independent entities: The first, which derive from uncharged spin— 1/2 electrons, can travel at 
different velocities when compared with the second, built from charged spinless electrons. Predicted 
theoretically in the early sixties, the spin-charge separation has attracted renewed attention since 
the first evidences of experimental observation, with usual mentions as a possible explanation 
for high-temperature superconductivity. In one-dimensional (ID) model systems, the spin-charge 
separation leads the frequencies of Friedel oscillations to suffer a 2k f — > Akp crossover, mainly 
when dealing with strong correlations, where they are referred to as Wigner crystal oscillations. 
In non-magnetized systems, the current density functionals which are applied to the ID Hubbard 
model are not seen to reproduce this crossover, referring to a more fundamental question: Are 
the Wigner crystal oscillations in ID systems non-interacting v-representable? Or, is there a 
spin-independent Kohn-Sham potential which is able to yield spin-charge separation? Finding 
an appropriate answer to both questions is our main task here. By means of exact and DMRG 
solutions, as well as, a new approach of exchange-correlation potential, we show the answer to 
be positive. Specifically, the v-representable Akp oscillations emerge from attractive interactions 
mediated by positively charged spinless holes - the holons - as an additional contribution to the 
repulsive on-site Hubbard interaction. 

PACS numbers: 71.10.Fd, 71.27.+a, 71.15.Mb 
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I. INTRODUCTION AND BACKGROUND 

Friedel oscillations are oscillations of the electron density that appear around inhomo- 
geneities and were first described by J. Friedel.- With renewed attention due to new possi- 
bilities of experimental realization, Friedel oscillations have emerged as a useful laboratory 
in the study of many-body systems,- 1 ^ in particular those described by the one-dimensional 
(ID) Hubbard chains. Contrary to higher dimensional systems, which use to be described 
by the Fermi-liquid theory, the ID metals assume non-Fermi-liquid behavior, and belong to 
a special class of Tomonaga-Luttinger liquids. 4,5 

In one dimension, and in second-quantized notation, the Hubbard model 6 (HM) is defined 

as 



where L is the number of sites, t the amplitude for hopping between neighboring sites, and 
U the local (on-site) interaction acting on site j. Occupation of each site is limited to two 
particles, necessarily of opposite spin. 

The density profiles of open Hubbard chains display Friedel-like oscillations, as it can be 
seen in Fig. [TJ obtained by means of Lanczos exact diagonalization and DMRG calculations. 
We note a clear change in the frequencies of oscillation upon variation of the local Hubbard 
interaction U. The increase of U is accompanied by an increase in the frequency, whose value 
is known to pass from 2k p at U/t = to 4kp in the limit of U/t — > oo, corresponding to the 
transition from three to six positive peaks in Fig. [T](a) [or from eight to sixteen in Fig. QJb)], 
where 2kp = tt(N + 1)/(L + 1), N is the number of electrons and hp is the Fermi wave 
vector Specifically, each positive peak in the density profiles is splitted into two new ones. 
The resulting 4k p Friedel oscillations, usually referred to as Wigner crystal oscillations,-" 1 ^ 
are connected with the spin-charge separation effect, where the one-dimensional confinement 
yields strongly interacting electrons which can break their spin and charge into two separated 
quasiparticles, the first built from uncharged spin— 1/2 and the second from charged spinless 
electrons.—"^ This mechanism, with recent evidences of experimental observation,—"— is 
usually mentioned as a possible explanation for high-temperature superconductivity.— ~— 

In the limit of U/t = 0, the density profiles n(j) of Friedel oscillations are exactly de- 
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FIG. 1: (Color online) Panel (a): Density profiles from Lanczos exact diagonalization, with 
L = 15 sites and N = 6 electrons. Panel (b): the same as before, but from DMRG calculations 
with L = 100 sites and iV = 16 electrons. The average densities were added to constant values in 
order to clarify the visualization. 
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whereas, in the limit of U/t — > oo, the Wigner crystal profiles are given by 

W (L + l) 2(L + l)sin(^)' {6) 

which is equivalent to the U/t = limit of spinless electrons. 

There have been many attempts to describe the 2k p — > 4k p crossover by means of a 
local-density approximation (LDA), or its spin-dependent version (LSDA), which includes 
to consider magnetized situations with ^ A^, 19 and the break of spin symmetry by 
means of a vanishing magnetic field when iV-f = N^M ~— Generalized gradient approximations 
(GGAs)^ and self-interaction corrections (SICs)^ have also been considered, with no success 
for non-magnetized systems. The question we intend to discuss here is: Is it possible, 
in a non-magnetized system, to reproduce all ranges of frequencies via a standard Kohn- 
Sham (KS) density-functional theory (DFT) calculation? In other words, we face the v- 
represent ability question: Is it possible that the strongly interacting 4kp regime - and the 
spin-charge separation - cannot be yielded by a spin-independent KS potential? 
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A. The concept of v-representability 



The Hohenberg-Kohn theorems 2 ^ state that the total energy of a many-body system is a 
functional of density, usually written as 

E[n] = T[n] + V ee [n] + V[n], (4) 

standing for the kinetic, interaction and external potential terms, respectively. The exact 
expression for the ground-state kinetic energy can be written as 25 

T M = E E f k '° (Vvl*iVv)> ( 5 ) 

c=t4 fc 

where t, ipk,a and fk, a are the kinetic energy operator, spin orbitals and their occupation 
numbers, respectively. The Pauli exclusion principle requires that < fk >(7 < 1. The total 
electron density is given by 

n ^= EE /j^IiM*)! 2 , (6) 

<T=t>4- k 

under the constraint of ^<x=t i ^fe " f k > a = N-^ + = N . For interacting systems, with 
fractional occupation numbers for /fe )(J , there are an infinite number of possibilities in Eqs. 
(jHJ) and ([6]). Kohn and Sham considered /fe i<T = 1 for the first N a orbitals and fk,a = for the 
rest, that is, they approximated T[n] by the non-interacting kinetic energy, 7s [n], exactly 
described by these choices of occupation numbers. In practice, V ee [n] is also approximated, 
by the classical Hartree expression, Esin], and all corrections to both approximations are 
considered in the exchange-correlation (XC) term, E^n]. 

The use of 7g[n], however, leads to a restriction on the density: It needs to be non- 
interacting v-representable, that is, there must exist a non-interacting ground-state with 
the density n(r). In other words, there must exist a XC potential Uj; C [n](r) = SE xc [n]/Sn(r) 
whose ground-state, in a non-interacting KS calculation, leads to the correct interacting den- 
sity profiles. There are examples in the literature of densities which are not non-interacting 
v-representable,— and a common procedure to circumvent this is to allow fractional occupa- 
tion numbers for f^^. The v-representability is not a conceptual problem to the Hohenberg- 
Kohn theorems, since it comes from the choice of Ts[n] for the KS approach of DFT. We 
can mention alternative approaches, as the orbital-free DFT, 27 which aims to approximate 
T[n] as an explicit functional of the density, and does not make use of the KS equations. 
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II. RESULTS 



In order to prove the v-representability, we can (a) use accurate solutions of ID Hubbard 
chains, such as the DMRG, and invert the KS equation to find the accurate XC potential^ 
which yields the accurate densities; or (b) find an XC potential, or propose a new one, which 
yields the correct oscillations. In either case, if the potential exists, this will be a numerical 
(and sufficient) proof that the 2kp —> 4kp crossover is non-interacting v-representable. We 
here will go through these two options. 

A. Accurate exchange- correlation potentials 

Using accurate density profiles as input, we numerically traced back the KS and XC 
potentials. The inversion of the KS equations is performed in a self-consistent cycle, as 
follows: Step 1: Introduction of a initial guess for the entire KS potential V^ s (j), for example, 
which contains only the Hartree term of a uniform system. Step 2: Determine the density 
profiles n°(j) obtained from V^ s (j). StepS: Compare n°(j) with the accurate density profiles 
n Accuratc (j). Step 4: If n°{j) < n Accurate (j), V£ s {j) = V£ s (j) - 5; if n°{j) > n Accurate (j), 
^Ks(i) = ^ksO) + with 5 < 1 x 10~ 5 . The 2 — > 4 cycle is repeated 7 times until 
n f (j) « n Accurate (j) ; with an accuracy of ~ 0.1%. Results are shown in Fig. [21 where we 
plot the accurate XC and KS potentials resulting from the density profiles of Fig. [T] and 
from additional data. 

The answer about the v-representality is positive, that is, it is possible to find spin- 
independent XC potentials which yield the accurate density profiles and their crossover in 
frequency. The same holds for the entire KS potential, fKs(j) = VH[n](j)+v xc [n)(j)+v ex t(j), 
with VH[n](j) = Un(j)/2 and here v ext (j) = 0. The comparison between Figs. [T] and [^a)-(d) 
shows that, even though a change in frequency is present for the densities, the XC potentials 
always oscillate with the same frequency, 4k f (identified by the six or sixteen negative peaks, 
for each number of electrons). The only changes are in the oscillation amplitudes, which 
are increased as the interaction U goes in the same direction. This behavior is specially 
unexpected for weak interactions, where the density profiles oscillate with 2k p, whereas 
the XC potentials do not follow the same path. As consequence of an increment in the 
XC amplitude, the entire KS potential suffers a transition from 2kp to 4k F , driving the 
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FIG. 2: (Color online) Panels (a)-(b) and (e)-(f): Accurate (on Lanczos exact diagonalization 
accuracy) XC and KS potentials obtained by means of the KS equation inversion: L = 15 sites and 
N = 6 electrons, with N+ = Ni. Panels (c) and (d): the same as before, but on DMRG accuracy 
for L = 100 sites and N = 16 electrons. The average potentials were added to constant values in 
order to clarify the visualization. 
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same path to the density, as seen from Fig. HJ Exceptions are some cases of very weak 
interactions, as pointed out in Figs. EJje) and (f), where the XC potentials may oscillate 
with an intermediate frequency 2k F < f < 4k F , identified by five negative peaks, and with 
no effects to the densities (which oscillate with 2kp for these cases, due to the very reduced 
XC amplitudes). 

As a first conclusion, the key ingredient which yields the 2kp — > 4k,F crossover in the 
density profiles is the change in amplitude of a 4kp oscillating XC potential. With the 
prominence of correlation effects, the spin-charge separation emerges from the XC potential, 
which oscillates with the same frequency as spinless electrons do. 

B. Approximated exchange-correlation potentials 

Once it is v-representable, the next option concerns the possibility of finding an approx- 
imated XC potential which yields the 2kp — > 4kp crossover. The basic Hohenberg-Kohn 
and Kohn-Sham theorems of DFT also hold for the Hubbard model, once the density n(r) 
is replaced by the on-site occupation number- 



In terms of this variable, local- (spin)-density approximations for Hubbard chains and rings 
have been constructed. Successfull examples are the analytical parametrization proposed by 
Lima et al.,— which is now called BALDA/LSOC, or its alternative route, which uses a direct 
numerical solution of the Bethe Ansatz integral equations, and is known as BALDA/FN 
(fully numerical) .— Both approaches can be employed in the usual way, via a KS calculation. 

The LDA functional, with n^(j) = n^{j), is not seen to describe the Friedel oscillations 
correctly, at least in the limit of Wigner crystallization. We here will not consider the 
application of any vanishing magnetic field, which can break the spin symmetry and yield 
the 4kp oscillations, however, at expense of a magnetized final solution with n^(J) 7^ n±(J) for 
iV-f = N±. Instead, we propose a spin- independent XC potential, based on a complement of 
the previous BALDA approaches. The total Hamiltonian of an one-dimensional interacting 
system is known to separate into two independent terms, of spin and charge, as follows: 11 




(7) 



<7=t,l 



H = Hq + Hp + H p = H + Hi, 



(8) 
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FIG. 3: (Color online) Schematic occupied states of three systems: (a) Non-interacting, with spins 
and charges together, (b) Strongly-interacting, with spinons and chargons following the fractional 
statistics of semions. (c) Non-interacting Kohn-Sham system, with spins and charges together 
under influence of holons (the chargon antiparticles) . 

where Hq, Hp and H p stand for the kinetic, spin and charge terms, respectively. Here, 
the spin densities (3(j) are built from uncharged spin— 1/2 electrons (spinons), whereas the 
charge densities p(j) are built from charged spinless electrons (chargons). Spinons and 
chargons are semions, that is, particles with statistics half that of regular fermions,— ~— 
with occupied states following the schematic representation of Fig. |3](b), where charge and 
spin were fractionalized to form the spin— 1/2 spinons (as a sum of two 1/4 spins) and 
the spinless chargons (as a sum of two — e/2 charges). In this context, considering spinons 
and chargons as distinct separated entities, the total interacting density can be written as 
n U) = /5(j)+p(j), with spin and charge allowed to travel at different velocities. Particularly, 
in the strongly- interacting limit, it has been shown that in some cases the spin densities can 
be considered as almost static in comparison with the charge moviment. 35 

In contrast to spin— 1/2 electrons and spin— 1/2 holes in the Fermi-liquid theory, excita- 
tions in one-dimensional systems are described in terms of spinons and positively charged 
spinless holes - the holons (the chargon antiparticles). It has been argued that spinons and 
holons attract each other at short distances.— In this context, as shown in Fig. E^c), we 
here propose that the occupied states of a non-interacting KS system are built by retaining 
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spin and charge together, at expense of the presence of holons, whose densities are given by 
p + (j). The KS potential can thus be written as a functional which depends on n KS (j) and 



with n KS (j) = n(j) and 

N 1 

p + ^)= EE o \M)\ 2 ' ( 10 ) 



CT=t,4 

that is, the holon density p + (j) is constructed by using the first N a unoccupied KS orbitals. 
This conjecture is specially conceived to deal with strong interactions, where the presence 
of holons will be felt by the electronic density as an attractive external potential. Based on 
Eqs. © and (ITU1) . we introduce the present approach: 

^Ks sent N(j) = v ext (j)+v H [n)(j)+v:r o3: [n)(j) 

-v H [ P + ](j) - 

= v ext (j) + v H [n](j) + v p x r n >]U), (11) 

with 

at ct 

»o*)= EE /^i^0')i 2 , (12) 

and the KS choices for j^ a . As approximated functional (approx) we chose the BALDA/FN. 
We include the holon density only in the XC potential, and use it in the KS calculation. A 
welcome feature of the present approach is the absence of the one-electron self-interaction 
error: In a non-magnetized LDA calculation, VfetO") = ^hU)- Therefore, for one-electron 
systems, n(j) = p + (j) = \i^ka{j)\ 2 , and then the present KS potential is exact, correcting 
the spurious self-interaction of one electron with itself. The results for BALDA/FN and 
the present correction are shown in Figs. H] and [5J for the electronic density n(j) and XC 
potential. Remember: BALDA/FN takes into account only the first two terms of Eq. ([9]). 

From Fig. H(a), while BALDA/FN yields the correct behavior for weak interaction, the 
2k p — > ^kp crossover is not recovered as U is increased, following the same trends already 
observed in the literature.— - — The same occurs with the BALDA/FN XC potentials of 
Fig. H](b), which oscillate in 2kp for all values of U (identified by the eight negative peaks). 
From Fig. H](c), the densities from the present approach also oscillate correctly for weak 
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FIG. 4: (Color online) Panels (a) and (b): BALDA/FN density and XC potential profiles. Panels 
(c) and (d): Present approach profiles. L = 100 sites and N = 16 electrons, with = N±. 

interactions, but, different from BALDA/FN, we do observe the 2kp — y 4kp crossover as U 
is increased, even though the amplitudes are subtly higher than the DMRG prediction of 
Fig. [TJ In addition, from Fig. 11(d), the XC potentials from the present approach oscillate 
in 4kp for all values of U, the only changes are in the amplitudes, the same as observed for 
the accurate XC profiles of Fig. |2j The spinless term of Eq. fl9]) was able to incorporate the 
key ingredient we mentioned before, i.e., the spin-charge separation by means of a change 
in amplitude of a 4k p oscillating XC potential. 

It is known that the prominence of correlation effects is not only driven by the interaction 
U, but also by low densities. For this reason, we have also considered Hubbard chains with 
different average densities n = N/L. Specifically, in Figs. E^a)-(b) we plot the LDA XC 
potentials for two values of n. Interestingly, for each case, there are values of U for which 
the LDA XC potentials oscillate in 4kp, as signed by the four (when = 4) or six (when 
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N = 6) negative peaks. As U is increased, the LDA XC potential suffers an incorrect 
transition to 2k f (identified by two or three negative peaks), and then yields the incorrect 
density profiles for larger U (not shown). This feature can be better understood by means 
of the XC potentials of homogenous extended systems, which are used in the construction 
of the LDA XC functionals, and are plotted in Figs. E^a)-(c). We note a common behavior: 
Near the average densities of the systems of Fig. [5j the homogeneous XC potentials undergo 
minimum values, whose positions tend to move to the right as the interaction is increased. 
For a given value of U, the LDA XC potential emerges to oscillate in 4kp only for the 
average densities surrounding these minimum values. Without the inversion points, as seen 
for U > 3t in Fig. E(c), an increase of n(j) will be always followed by a decrease of v xc [n](j), 
and vice- versa, explaining correspondence between Figs. HKa) and (b). In Fig. E(d) we 
plot the values of n < 1 for which the LDA XC potentials oscillate with 4/sp: For average 
densities above this curve, the LDA XC potentials oscillate with an intermediate frequency 
2/cp < f < 4kp, as identified by three and four negative peaks in Figs. ISfa) and (b), 
respectively, in accordance with the similar behavior of Fig. E|e). Bellow the curve, the 
LDA XC potentials oscillate with the wrong frequency 2fcp. 

To conclude, (i) it should be stressed that even though the LDA XC potential may 
oscillate with / and 4k f, these frequencies occur for values of U in which the XC amplitudes 
are so reduced that the density profiles are not affected, retaining the 2kp oscillations in 
all cases we have implemented; (ii) the curve of Fig. [6]^d) can be also employed to mark 
the necessity of using the KS states exposed in Fig. EJ^c), where the spin-charge separation 
emerges from the inclusion of spinless holons. As expected, the separation is prominent 
both for strong interactions and/or low densities, situation indicated in Fig. [6]^d) as the 
"4fcp region". 

The present approach yields the correct frequencies of density profiles for all values of U 
and n we have investigated, including the systems of Fig. [5] (not shown). This is another 
numerical proof that is possible to find a spin-independent KS potential which yields the 
Friedel and Wigner crystal oscillations, for all ranges of interaction and prominence of corre- 
lation effects. All this in the following implementation scenario (on a supercomputer): The 
computing time of DMRG scales with ~ 10 1 hours, while the present DFT approach with 
~ 10° minutes. We have also considered other flavors of self-interaction corrections, such as 
the full Optimized Effective Potential (OEP) implementation of the usual Lundin-Eriksson 37 
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FIG. 5: (Color online) BALDA/FN XC potential profiles: (a) L = 17 sites and iV = 4 electrons; 
(b) L = 13 and N = 6. Both cases with iV t = N±. 




FIG. 6: (Color online) Panes (a)-(c): XC potential profiles for extended homogeneous systems 
obtained by means of the Bethe Ansatz integral equations. Panel (d): Values of n for which the 
homogeneous XC potential undergoes a minimum - and the LDA XC potential oscillates with 4fcp 
(for n < 1). 



12 



and Perdew-Zunger— approaches. However, none of them were able to yield the Akp oscil- 
lations, what teach us that the 2kp — > 4kp crossover, or the spin-charge separation, is not 
only a problem of one-electron self-interaction error.— More than correcting the spurious 
self-interaction of LDA, which tends to delocalize electrons, an accurate functional should 
include the physics of positively charged holons. 



C. Analytical conjecture 

Based on Eqs. and ([2]), it is possible to propose an analytical conjecture to describe 
the Friedel and Wigner crystal oscillations, valid for any value of U. It is usually written 
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as: 



"W-(TTT)-^ [3in(A)] <^-^ [sin(A)] - ■ (13) 

The parameters A\, A 2 and K c depend on U and are related to the amplitudes and oscilla- 
tions decay. Specifically, Ax and A 2 are usually referred to as Friedel and Wigner amplitudes, 
respectively. K c is the correlation exponent, one of the Luttinger parameters.— A similar 
expression has been proposed for the density distribution in a harmonic trapping potential.— 
We here fitted the density profiles of Figs. d](b) and H] to the expression f[T3"j) . with the values 
of Ax, A 2 and K c expressed in Table Has a function of U. In comparison with DMRG, LDA 
fails to predict the correct values of A 2 , as well as, the correct decay of A\ and K c . The 
present approach performs much better, mainly for the K c values, which are obtained in 
great agreement with DMRG, even though it is seen to overestimate the weight of A 2 for 
large U, as expected from the density profiles of Fig. HKc). 

Beyond the local interaction U, all parameters in Table [I] depend on the number of 
particles N, chain size L and average density N/L. For this reason, these values cannot 
be considered as universal - they are only a picture of a specific system with N = 16 and 
L = 100. A more complete analysis follows in Fig. [7J where we plot the crossover points 
with A\ = A 2 in the U - N/L plane. Based on DMRG calculations, an accurate prediction 
for L = 200 is N/L 0.034 U. 7 - For a small chain of L = 15 sites, we here introduce an 
additional accurate expression, N/L w 0.058 U, based on Lanczos exact diagonalization. 
For average densities above these curves, systems display Friedel-like oscillations with 2kp 
{Ax > A 2 ). Below it, the Wigner crystal behavior with 4kp is prominent {Ax < A 2 ). The 
constant BALDA/FN curve at zero means that it predicts Ax > A 2 for all values of U and 
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TABLE I: Parameters A±, A 2 and K c obtained by means of a fit of expression (|13p to the data 
of Figs. Q] (b) and HI The nomenclature "x.xx / y.yy / z.zz" indicates the parameters values for 
U = 2t, U = At and U = \2t, respectively. 
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FIG. 7: (Color online) Friedel and Wigner crystal amplitudes. Crossover points with Ai = A 2 : 
For average densities above the curves, systems display Friedel-like oscillations. Bellow them, the 
Wigner crystal behavior is prominent. 

N/L. The curves from the present approach emerge to be much more accurate, mainly for 
low densities and/or moderate interactions. Its great improvement over BALDA/FN is due 
to a more accurate description of the 2kp 4:kp crossover. For strong interactions and 
higher average densities, as seen for L = 200, the present approach tends to fit below the 
accurate prediction, yielding a larger Friedel-like region. 
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III. CONCLUSIONS 



We have dealt with the v-representability issue of Wigner crystal oscillations in one- 
dimensional Hubbard chains. Initially, we posed two questions concerning non-magnetized 
systems: (1) Are the ID Wigner crystal oscillations non-interacting v-representable? (2) Or, 
is there a spin-independent KS potential which is able to yield the spin-charge separation? 
We have concluded that the answers to both questions are positive, that is, we do not face 
a v-representability problem. Different from some previous studies, we have not considered 
the break of spin symmetry, which is known to yield the 2k p — > 4kp crossover. Instead, 
by means of accurate many-body solutions for small and extended chains, we inverted the 
KS equation, showing that it is possible to find a spin-independent XC potential which 
reproduces the accurate densities. 

In addition, we have proposed a new XC potential, which incorporates, by means of 
unoccupied KS orbitals, the physics of positively charged spinless holes (holons). Specifically, 
holons act on mediating attractive interactions, as an additional contribution to the repulsive 
on-site U. The present approach is also able to reproduce the 2kp — > Akp crossover, either 
for small or extended chains, with correlation exponents obtained in great agreement with 
DMRG data. 

It remains to be seen whether the present spin-independent XC potential ( TXT]) can cir- 
cumvent known failures of LDA,— as well as, be successfully applied to other classes of 
systems, with impurities and external confinements, as the example of ultracold fermions in 
one-dimensional traps.— 
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